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PACS 64 . 70 . D- - Solid- liquid transitions 

PACS 6 1.50. Ah - Theory of crystal structure, crystal symmetry; calculations and modeling 

Abstract. - Phase field crystals (PFC) are a tool for simulating materials at the atomic level. 
They combine the small length-scale resolution of molecular dynamics (MD) with the ability 
to simulate dynamics on mesoscopic time scales. We show how PFC can be interpreted as the 
result of applying coarse-graining in time to the microscopic density field of molecular dynamics 
simulations. We take the form of the free energy for the phase field from the classical density 
functional theory of inhomogeneous liquids and then choose coefficients to match the structure 
factor of the time coarse-grained microscopic density field. As an example, we show how to 
construct a PFC free energy for Weber and Stillinger's two-dimensional square crystal potential 
which models a system of proteins suspended in a membrane. 



Introduction. — Phase field crystals (PFC) were in- 
troduced in [1] as a tool for performing long-time simula- 
tion of materials at the microscopic scale. The approach is 
an example of the phase field approach to non-equilibrium 
modeling and simulation. The system of interest is rep- 
resented by a field corresponding to the atomic density of 
the material. Its dynamics are simulated dissipatively ac- 
cording to free energy minimization with additional noise. 
The advantage of PFC over other phase field models is 
its ability to naturally resolve microscopic structure: the 
PFC field corresponding to a solid in equilibrium has a 
periodic structure matching that of the crystal stucture 
of the material. The advantage of PFC over molecular 
dynamics (MD) is that fast degrees of freedom are not 
present in the model. This feature allows efficient simula- 
tion over mesocopic time scales. The PFC approach can 
model elastic and plastic deformation of crystals as well as 
the liquid-solid transition [2] and diffusion of defects [3]. 

Previous work [2] demonstrated how phase field crys- 
tals could be constructed that matched the first peak in 
the structure factor of a two-dimensional hexagonal crys- 
tal material and demonstrated their use for the simulation 
of a variety of non-equilibrium phenomenon. Here we will 
accomplish two things: (i) show that phase-field crystals 
can be interpreted as the result of applying coarse-graining 
in time to microscopic molecular dynamics (ii) show how 



to generalize phase-field crystals in order to model dif- 
ferent materials and match more than the first peak in 
the structure factor. The system we will use as an exam- 
ple for this is Weber and Stillinger's [4] two-dimensional 
square-crystal forming system. This system was originally 
proposed as a model of various systems where square two- 
dimensional crystals occur, including proteins suspended 
in a membrane and planar levels in colloidal suspensions. 
We will show that we can obtain a PFC free energy that 
matches arbitrarily many peaks of the stucture factor of 
the time-coarsened microscopic density for this potential. 

Coarse-Graining Molecular Dynamics. — We 

consider a system of particles in the plane each with 
mass m and with positions and momenta Pi^ i = 
1, . . . , A^. The energy of the system is given by 
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(1) 



The V used by Weber and Stillinger [4] contains both two- 
body and three-body terms: 



V(ri, . . . ,rAr) = 



(2)/ 



,-|) + A^^(^)(r,,r„rfc). 
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The pair potential is given by 

exp[(r - a)~^], < r < a 



A{r- 
0, 



-12 



r > a, 



where A = 6.767441, a = 2.46491832. The three-body 
term is given by 



where rij = \ri — rj\ and 6i is the angle between r j — 
and r/c — r^. The function is 

h{r, s, 0) = exp[(r - aa)"^ + - ag)"^] sin2(2l9), 

when r, s < as, and is zero otherwise. Weber and Still- 
inger choose as = 1.8 and A = 5. Without the three- 
body terms the potential is qualitatively the same as the 
Lennard- Jones potential, and the system has a hexagonal 
crystal phase like that modeled by the PFC in [2,3]. The 
addition of the three-body term makes hexagonal struc- 
ture unfavourable and the crystal phase has square sym- 
metry [4]. The dynamics of the system were simulated 
using molecular dynamics. 

The microscopic density field is given by pmir^t) = 
X]i=i ^(^ — ri{t))^ where ri{t) is the position of the zth 
particle at time t. We consider the density coarse-grained 
in time: ^ 

p{r) = - I dt p^(r) (2) 

If the original system is in the liquid phase, for large r 
the field p will be equal to a constant plus small, spatially 
homogenous fluctuations. If the original system is in the 
solid phase, p will have a steady periodic form for large 
r corresponding to the structure of the crystal. Fig. 1 
shows density plots of p of a system in solid phase for 
four different coarse-graining times r. As r is increased 
fluctuations about equilibrium become smaller and occur 
over a smaller time scale. 

The Phase Field Crystal. — Our goal is to obtain 
a phase field model for the dynamics of the coarse-grained 
density field p given by eq. (2) for large but finite r. This 
goal can be broken into two distinct parts. Firstly, we need 
to determine a free energy functional whose minimum is 
a good approximation to the equilibrium configuration of 
p. Secondly, we need to specify evolution equations that 
describe how the system approaches equilibrium. 

Free Energy. The motivation for the particular free 
energy functional T we use comes from classical den- 
sity functional theory [5,6]. The goal of Ramakrishnan 
and Yussouff 's density functional theory is to estimate a 
free energy for which the infinite-time averaged den- 
sity p = poo = {pm) is a minimum. They derive this free 
energy formally by integrating the partition function; the 
calculation can only be completed in the case of the ideal 
gas. To obtain usable approximations they postulate that 
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Fig. 1: Grey-scale plot of time- averaged density of portion of 
Weber-Stillinger MD system for reduced temperature 0.6, re- 
duced density 0.77277. The time intervals for averaging are 1, 
10, 100, 1000 time units. 

p has a free energy J^{p) = J^id{p) + J^ex{p) where J^id is 
the ideal gas free energy and J^ex is a functional Taylor 
series expansion in p about the uniform liquid state den- 
sity [5]. Typically, only two terms in the expansion are 
taken for most applications [7]. 

Our goal is somewhat different than that of density func- 
tional theory; we are interested in densities averaged over 
finite times r, rather than infinite times. This allows us to 
model non-equilibrium phenomena. We start with free en- 
ergies of density functional theory form and then add noise 
to account for our merely finite length of time-coarsening 
interval r. Moreover, rather than deriving coefficients in 
the free energy from liquid state equilibrium correlation 
functions, we choose coefficients that will be easy to com- 
pute, give the correct qualitative behaviour and have pa- 
rameters that can then be tuned to match MD. 

Our free energy is of the form 



Hp) 



drT p{r){\og{p) — 1) — / drp{r)p 



2 



(iri(ir2 C2 (ri , r2 )p(ri )p(r2 ) 

+ T" /// ^^i^^2C^^3C3(ri,r2,r3)p(ri)p(r2)p(r3) 

The last cubic term is necessary to obtain a square crys- 
tal structure and is analogous to the three-body potential 
term in [4]. The variable T corresponds to temperature. 

According to density functional theory, the kernal C2 
should be a multiple of the direct pair correlation func- 
tion of the system in the liquid state [8] . Following this we 
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choose our C2 to be radially symmetic and translationally 
invariant, giving C2(ri,r2) = c(|ri — r2|) for some func- 
tion c that decays to zero at infinity. However, instead of 
matching c to the correlation function exactly we instead 
choose a simple, compact bump function. Our results are 
quite insensitive to the exact details of the choice of c we 
use, the most important quantities being the height and 
the location of the peak of the Fourier transform of c, both 
of which can be controlled by scaling. This is predicted 
by the early papers on density functional theory [5] and in 
later analyses tying density function theory to phase field 
crystals [9]. The function c(r) = 1 for |r| < ro, c(r) = 
otherwise works for our purposes, but in order to minimize 
grid artefacts we use the following smoother choice: 

If arctan 20(r/ro — 1) 
"^'^ = 2[ 2^ ' 

For ease of computation, this function is then truncated 
at radius 2ro and shifted down to be continuous. We fix 
ro = 6 and /c^i = 1. Again, the C3 from density functional 
theory is a three-point correlation function that could in 
principle be computed from molecular dynamics. Instead 
we choose a simple, smooth, compactly supported C3 that 
gives the qualitatively correct behaviour: 

C'3(ri,r2,r3) = U {r 12, r 13,0 123) ^ U{r 21, r 23,0 213) 

+ ^(^31,^32,6>3i2) 

where Vij = \rj — and Oijk is the angle between rj — Vi 
and r/e — r^. We let 

t/(ri2,r23,^) = i/(ri2)i/(r23)sin'(2^) 



This yields dynamics in which mass is conserved. We will 
also consider 



dpir) 



-Tp 



•ra2 + ^2^^^(r,t) (4) 



dt ^ Sp{r) 

which corresponds to constant chemical potential. Here 
(C(r,t)C(r',t')) = S{r-r')S{t-t'). 

In both cases the noise should be interpreted in the Ito 
rather than the Stratonovich sense. The Fa^ term in 
eq. (4) does not represent an external force but is a con- 
sequence of using multiplicative noise. It is necessary for 
the dynamics to satisfy detailed balance, as can be checked 
with the Fokker-Planck equation. Both of these dynam- 
ics yield an equilibrium distribution for p proportional to 

Naturally, neither of these equations are expected to 
correctly capture dynamics far from equilibrium and close 
to equilibrium they will only be accurate in situations 
where molecular dynamics is well matched by Brownian 
dynamics. The non-conserved dynamics has the same 
equilibrium properites as the conserved dynamics and can 
be simulated numerically more efficiently, so we shall use 
them for our studies here. We will show that for suitably 
chosen a and parameters in the dynamics of eq. (4) has 
an equilibrium distribution that matches that of the time- 
averaged microscopic density p, at a given temperature, 
density, and coarse-graining time r. 

The resulting evolution equations when we use non- 
conserved dynamics eq. (4) with given by (3) are 



where 



Hir) 



exp[(r 



dp 
dt 



-Fp [T log p - /i + kACp + ksVip, p)] 



+Fcr^ 



for Try 



< r < Tr, 



= 6, r„ 



and H{r) = otherwise. We fix ^j^^re Cp(ri) = J dr2c{ri - r2)p{r2) and 



= 10 and ks = 0.05. We did not find 



that the resulting equilibrium configuration was sensitive 
to the details of C3. The purpose of the third-order term is 
only to make a hexagonal pattern unstable and the square 
pattern more energentically favourable. 

Equations for time evolution. Now that we have a 
form for the free energy JF, we need to state equations for 
the time evolution of the phase field. In the present work 
we do not match non-equilibrium properties of the PFC 
to those of MD. Instead we only want plausible dynamics 
with the correct equilibrium properties. In this way we 
can use the dynamics to sample efficiently from the equib- 
librium states. Accordingly, we use a form of equation 
from the dynamic density functional theory of interacting 
Brownian particles [10]: 



dp{r) 
dt 



: rV • LoV 



5J^ 

5p{r) 



'■(^V2^^ff{r,t)'), (3) 



where 77 (r, t) is a mean zero noise field with 



{r],{r,t)r]j{r\t^)) 



(j^8ij8(r-r')8{t-t'). 



Here ^ is noise white in space and time. 

Numerical Results. — To numerically simulate the 
dynamics p is discretized on a rectangular grid. In the 
units we use, the distance between two adjacent points on 
the grid is 1 unit. All integrals are evaluated by straight- 
forward summation on the grid. The equation are inte- 
grated in time using the explicit Euler method, typically 
with a time step of length At = 0.01. 

For particular values of the parameter p and without 
noise (<j = 0) the field p has two possible equilibria selected 
depending on T. For large T the log term dominates and 
the stable equilibrium is a constant p corresponding to a 
liquid state. As T is decreased, this state becomes unsta- 
ble and the system has a periodic free energy minimum 
resembling the fourth plot in fig. 1. With the parameters 
described here the distance between two adjacent bumps 
in the periodic phase is approximately 8.5 units. 
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Fig. 2: Grey-scale plot of p on a portion of the computational 
domain during simulation with non-conserved dynamics. Plots 
shown at 0, 33, 66, 99 time units. 



Nonequilibrium simulation. Fig. 2 shows a non- 
equilibrium simulation of the system with non-conserved 
dynamics. With /i = 10, T = 2.2, a = the system 
was started in the equilibrium crystal state except for one 
one patch consisting of square of 25 particle bumps which 
is rotated 45°. The system was integrated in time until 
close to the uniform solid configuration equilibrium. In 
grey-scale we show density plots of the system for four 
representative times. 

Structure Factor. We now compare results numeri- 
cally obtained using this free energy with those of time- 
averaged molecular dynamics. For a wide range of tem- 
peratures, densities, and time-coarsening times of the MD 
system, we can choose T, /i and a for the PFC to obtain 
good agreement, as we discuss further below. 

Fig. 3 shows the structure factor for a time-averaged MD 
simulation together with the structure factor for an equi- 
librium PFC configuration where parameters have been 
chosen to obtain the best fit. The MD structure factor 
was taken from a 400 particle system in a perfect square 
crystal of density 0.77277 (leading to the minimum en- 
ergy for the ground state) and reduced temperature 0.55. 
The PFC was obtained from a non-conserved dynamics 
run with /i = 10, T = 0.8, a = 0.02. The values for these 
parameters were chosen as follows. Since vertical and hor- 
izontal scales on the structure factor could be matched 
trivially be rescaling the dimensions of density and dis- 
tance, we only concerned ourselves with matching the rel- 
ative heights of the peaks and their relative spacing. With 
(7 = (no noise) /i and T were first chosen in order to ob- 




Fig. 3: Comparison of structure factor S{q) versus wavenumber 
q for MD (line and dots) and PFC (circles). For clarity only 
some of the data points are shown for PFC. The inset shows 
the same data on a logarithmic scale. 



tain a square crystal pattern. This occurs for a broad 
range of these parameters and immediately gives the cor- 
rect spacing of peaks in the structure factor. No further 
adjustments of /i were necessary. The value of T was then 
adjusted to match the ratio of the first two peaks of the 
structure factor to that of the MD results. As T goes to 
zero (a perfect crystal) the relative height of the second 
peak grows, and as T goes to the melting point the second 
peak shrinks. As soon as the second peak was mathced all 
the remaining peaks agreed to the degree shown in Fig. 3. 
At this point the structure factor of PFC consists purely 
of peaks and it was necessary to increase a until the back- 
ground rose to match that of the MD structure factor's 
background. 

Melting Transition. An open question in the theory 
of phase transitions is whether the melting transition form 
solid to liquid in two-dimensional systems occurs via a sin- 
gle phase transition or if there are two phase transitions 
with an intermediate hexatic phase [11]. Weber and Still- 
inger studied this question for their system using MD, and 
found no evidence for a two-stage melting process [4] . We 
perform an analogous experiment using our PFC model of 
the system. As in [4], we examine the liquid-solid transi- 
tion by starting the system in a perfect square crystal and 
increasing the temperature (in our case, the parameter T) 
slowly. We used the non-conserved dynamics with /i = 30, 
and a = 0,0.01,0.02,0.04. A phase transiton should ap- 
pear in a plot of mean density versus temperature as a 
discontinuity in the solution or its derivative. In order 
to highlight any such points of non-smoothness, we first 
smooth the mean density as a function of T and then con- 
sider its derivative with respect to T. In fig. 4 we show the 
negative derivative of mean density p with respect to T as 
a function of T. We observe only one peak in the deriva- 
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[11] David R. Nelson, Defects and Geometry in Con- 
densed Matter Physics (Cambridge University Press) 2002, 
sect. 2.4 




Fig. 4: Plot of -dp/dT as a function of T for the PFC system, 
where p is mean density and T is temperature. Results are 
shown for four different values of the noise strength a. Plots 
are offset vertically from zero for clarity. 



tive for each run of the simulation, suggesting a one-stage 
melting process. Our results are consistent with those of 
Weber and Stillinger. 
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